home *** CD-ROM | disk | FTP | other *** search
/ Magnum One / Magnum One (Mid-American Digital) (Disc Manufacturing).iso / d18 / nrpas13.arc / DFPMIN.PAS < prev    next >
Pascal/Delphi Source File  |  1991-05-01  |  2KB  |  79 lines
  1. PROCEDURE dfpmin(VAR p: glnarray; n: integer; ftol: real;
  2.        VAR iter: integer; VAR fret: real);
  3. (* Programs using routine DFPMIN must supply a
  4. FUNCTION fnc(p: glnarray):real; and a
  5. PROCEDURE dfnc(p: glnarray; VAR g: glnarray);
  6. which evaluate a function and its gradient. They must
  7. also define the types
  8. TYPE
  9.    glnarray = ARRAY [1..n] OF real;
  10.    glnbyn = ARRAY [1..n,1..n] OF real;
  11. in the main routine. *)
  12. LABEL 99;
  13. CONST
  14.    itmax=200;
  15.    eps=1.0e-10;
  16. VAR
  17.    j,i,its: integer;
  18.    fp,fae,fad,fac: real;
  19.    xi,g,dg: glnarray;
  20.    hdg: glnarray;
  21.    hessin: glnbyn;
  22. BEGIN
  23.    fp := fnc(p);
  24.    dfnc(p,g);
  25.    FOR i := 1 TO n DO BEGIN
  26.       FOR j := 1 TO n DO BEGIN
  27.          hessin[i,j] := 0.0
  28.       END;
  29.       hessin[i,i] := 1.0;
  30.       xi[i] := -g[i]
  31.    END;
  32.    FOR its := 1 TO itmax DO BEGIN
  33.       iter := its;
  34.       linmin(p,xi,n,fret);
  35.       IF ((2.0*abs(fret-fp)) <= (ftol*(abs(fret)+abs(fp)+eps)))
  36.          THEN GOTO 99;
  37.       fp := fret;
  38.       FOR i := 1 TO n DO BEGIN
  39.          dg[i] := g[i]
  40.       END;
  41.       fret := fnc(p);
  42.       dfnc(p,g);
  43.       FOR i := 1 TO n DO BEGIN
  44.          dg[i] := g[i]-dg[i]
  45.       END;
  46.       FOR i := 1 TO n DO BEGIN
  47.          hdg[i] := 0.0;
  48.          FOR j := 1 TO n DO BEGIN
  49.             hdg[i] := hdg[i]+hessin[i,j]*dg[j]
  50.          END
  51.       END;
  52.       fac := 0.0;
  53.       fae := 0.0;
  54.       FOR i := 1 TO n DO BEGIN
  55.          fac := fac+dg[i]*xi[i];
  56.          fae := fae+dg[i]*hdg[i]
  57.       END;
  58.       fac := 1.0/fac;
  59.       fad := 1.0/fae;
  60.       FOR i := 1 TO n DO BEGIN
  61.          dg[i] := fac*xi[i]-fad*hdg[i]
  62.       END;
  63.       FOR i := 1 TO n DO BEGIN
  64.          FOR j := 1 TO n DO BEGIN
  65.             hessin[i,j] := hessin[i,j]+fac*xi[i]*xi[j]
  66.             -fad*hdg[i]*hdg[j]+fae*dg[i]*dg[j];
  67.          END
  68.       END;
  69.       FOR i := 1 TO n DO BEGIN
  70.          xi[i] := 0.0;
  71.          FOR j := 1 TO n DO BEGIN
  72.             xi[i] := xi[i]-hessin[i,j]*g[j]
  73.          END
  74.       END
  75.    END;
  76.    writeln('pause in routine DFPMIN');
  77.    writeln('too many iterations'); readln;
  78. 99:   END;
  79.